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Abstract 

Tomography  is  well  known  because  of  its  many  applications.  Although  theoretically 
solved,  the  numerical  implementation  of  tomographic  reconstruction  algorithms  is  still  a 
difficult  problem.  In  this  article  the  numerical  implementation  of  a reconstruction  method 
using  Cesaro-means  and  Newman-Shapiro  operators  is  described.  The  key  point  herein 
is  the  use  of  suitable  quadrature  formulae  on  the  sphere.  It  turns  out  that  in  the  context 
described  product  Gaussian  formulae  are  best  suited.  The  algorithm  is  tested  at  the  so 
called  Shepp-Logan  phantom  which  is  a three  dimensional  model  of  a human  head. 


1 Introduction  and  notation 

The  problem  in  tomography  is  to  reconstruct  a function  F from  its  Radon  transform  suf- 
ficiently well.  Since  certain  classes  of  functions  can  be  expanded  into  series  of  orthogonal 
polynomials  it  is  essential  to  exploit  the  action  of  the  Radon  transform  on  orthogonal 
polynomials  and  on  polynomials  in  general. 

This  approach  is  the  more  interesting  since  the  inverse  of  the  Radon  transform  for 
polynomials  is  known  explicitly. 

The  convergence  of  orthogonal  expansions  to  the  given  function  is  often  achieved  only 
by  applying  a summability  method.  The  application  of  such  methods  can  be  interpreted 
as  a kind  of  “filter  technique”  which  is  necessary  for  sufficiently  good  reconstructions.  The 
combination  of  an  expansion  of  the  function  and  the  application  of  suitable  summability 
methods  leads  to  promising  reconstruction  algorithms. 

In  this  article  two  examples  for  a summability  method  and  their  implementation  are 
presented  — the  Cesaro-means  and  Newman-Shapiro-means.  After  some  introductory 
remarks  on  Laplace-series  at  the  end  of  this  section,  in  Section  2 the  theory  of  sum- 
mability methods  needed  here  is  presented.  In  Section  3 this  theory  is  applied  to  the 
reconstruction  of  functions  from  their  Radon  transform.  Section  4 describes  the  nu- 
merical implementation  of  the  reconstruction  formula  which  is  tested  on  the  so  called 
Shepp-Logan  phantom  of  a head  in  Section  5. 

In  this  article  the  following  notation  is  used.  Let  Br  denote  the  unit  ball  in  IRT , Sr~ 1 
denote  the  unit  sphere  and  Zr  :=  [-1, 1]  x 5r_1.  xy  denotes  the  Euclidean  product  of 
x,y  £ Rr. 
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The  spaces  of  restrictions  of  r-variate  polynomials,  homogeneous  polynomials  and 
homogeneous  harmonic  polynomials  of  degree  g £ INo  onto  a subspace  X c Mr  ( X = 

Sr~ 1 or  X = Br)  are  denoted  by  1P^(X),  IP £ (X),  H^  (X),  respectively.  The  space 
C{Sr~l)  of  all  continuously  differentiable  functions  is  provided  with  the  inner  product 
< F,G  >:=  fSr-i  F(x)G(x)dx.  The  surface  measure  of  the  sphere  is  denoted  by  ujt~i  —< 
1,1  >• 


Let  denote  the  Gegenbauer  polynomials  of  degree  g and  index  A and  C * = 
Ch/Cp(  1)  the  normalized  Gegenbauer  polynomials.  The  reproducing  kernel  function  of 
* 2 LL  -j-  T 2 r— 2 

JHT  (Sr~l)  is  given  by  G^_(xy)  = - — CM 2 (xy),  the  normalized  reproducing 


kernel  Gfl  is  defined  by  G,, 


(r  - 2)u>r- 

GJG,(  1). 


Let  Y £ {C(Sr~1),  L2(Sr~1),  Lp(Sr~1)}.  For  f EY  let 

oo  oo  „ 

L{f,x)  = Y (A vf){x)  - Y / f(v)Gu{xy)dy 
„=o  »=oJS’-' 


(1.1) 


be  the  Laplace-series  of  /,  where  (Avf)(x)  :=  Jsr_  1 f(y)Gu{x,y)dy  is  the  orthogonal 
* 

projection  of  / onto  HI  (Sr~1)  and  the  partial  sums  L^(f,  x)  = (A vf)  (a:)  are  the 

orthogonal  projections  of  / onto  WJl(Sr~1}. 

Whereas  for  Y = L2(Sr~1)  it  is  known  that  the  partial  sums  LJf,x)  converge  to  / 

T 

in  norm,  no  convergence  is  obtained  for  Y = C(ST~1)  or  Y — Lp(Sr~1)  for  p > 2H — 

r — 2 

2 

and  p < 2 (see  e.g.  [l]p.211).  Applying  a summability  method  the  situation  changes. 

r 

2 Summability  methods 

Let  A = (aMv)n^£iN0  be  an  infinite  matrix  for  which  the  elements  £ H fulfil  the 
following  properties. 

(i)  = 0 for  v > g, 

(ii)  lim^oo  a ^ = 1 for  v £ {0, 1}, 

(iii)  K m(£)  > 0 for  -1  < £ < 1,  where  K ^ :=  oan-vGv. 

If  with  the  aid  of  a summability  method  the  kernel  Gv  in  (1.1)  is  substituted  by  a kernel 

t1 

K^  = Ya^Gl'  (2-1) 

i/=0 

then  the  operator  LA  defined  by  the  transformed  series 

LA(f,x)  = lim  f f(y)Kl_l(x,y)dy  (2.2) 

can  be  shown  to  converge  pointwise  to  the  identity  provided  that  for  the  kernel  Kp  the 
properties  (i)— (iii)  of  the  matrix  A are  valid. 


U.  Maier 


Remark  2.1  The  coefficients  atll,  can  be  obtained  from 

<V  = (£^Gv(f.))(t)  = f Gv{tx)Ktl{tx)dw{x),  t € S' 
Jsr -1 


r — 1 


For  A being  the  matrix  of  the  Cesaro-means  the  proof  was  given  by  Kogbetlianz  [4] 
first.  Berens  et  al.  [1]  give  a proof  for  Cesaro-means  as  well  as  for  Abel-Poisson-means. 
They  also  prove  results  on  the  order  of  convergence  and  the  corresponding  saturation 
classes.  The  convergence  proof  for  Newman-Shapiro  operators  (V  = C(Sr~1))  can  be 
found  in  Reimer  [7]. 


2.1  Cesaro-means 


For  Cesaro-means  the  coefficients  afU,  in  the  summability  method  have  to  be  chosen  as 


(l)<i  (fc  + 

(*  + !)„  (l)M-„ 


(2.3) 


where  {jp)q  = p ■ (p+ 1)  • . . . • (p+ q — 1)  denotes  the  Pochammer  symbol.  Then  the  kernels 
Kfj,  in  (iii)  take  on  the  form 


K,= 


(!)m  (fc  + l)^ 


(2.4) 


Convergence  of  the  transformed  Laplace-series  (2.2)  is  valid  for  k > (r — 2)/2;  for  k > r— 1 
the  operators  even  are  positive  (see  Kogbetlianz  [4]). 


2.2  Newman-Shapiro  summability  method 

In  [8]  Reimer  considers  kernel  polynomials 


•K2t/+i(£)  :=  K2u(0  :=  pi^+i 


Gu+i(0 


as  used  by  Newman-Shapiro  [5].  Here,  pv+i  is  the  largest  root  of  Gv+\  and 


9u+ i = {r-  2)wr_i 


1 - Vl+i  fv  + r-  2 


(2u  + r)2  \ r — 3 


1 -vl 


+1 


1 


2v  + r Gi,-j-i(l) 


(2.5) 


(2.6) 


The  coefficients  a„„  in  the  Newman-Shapiro  operators  can  be  calculated  to  be 


v v 


9v+i 

j= o 1=0 


(2i/  + r)2 

(^  + 1)2 


Gj{Pi/+i)Gi(rju+i) 

(G,(^+1))2 


(j  + A)(/  + A) 

U)r-i\2 


* Wk  Wj-k  (A )i-k  (l)j+;-2fc  (2A )j+i~k  _ c 
fr'g  (!)fc  (!)i-fc  (l)i-fc  (2A)j+;_2fc  (A  + l)j+i-k  u'3+l  2k 


(2.7) 


where  6vj+i- u denotes  the  Kronecker  delta  and  A — 2-=^. 

The  matrix  A defined  by  the  Newman-Shapiro  operators  fulfils  the  properties  (i)-(iii) 
(see  Reimer  [8]). 
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Remark  2.2  The  corresponding  partial  sum  operators  are  nonnegative  with  positive 
a^.  For  continuous  and  differentiable  functions  even  more  is  valid  (see  Reimer  [ 8 ]): 
whereas  for  continuous  functions  the  approximation  error  is  of  order  Of/i”1)  , functions 
F G C:i(Sr~1),  j € {1,2},  have  an  error  of  order  0(p~j). 

3 Application  to  tomography 

The  Radon  transform  TZ  : C(Br)  — > C(Zr)  is  defined  by 

(KF)(s,t)~  J F(st  + v)dv,  FeC(Br),  ( s,t)  e Zr , (3.1) 

v.Lt 

V^<1  — 8^ 

which  means  that  the  Radon  transform  77.  of  F is  determined  by  integrating  F over  all 
hyperplaneS  of  dimension  r — 1.  This  map  can  also  be  defined  for  functions  in  I?  (JRr), 
L2(Br),  the  Schwartz  space  S(lRr ) or  some  Sobolev  spaces.  7 Z is  continuous  on  all  of 
these  spaces,  whereas  the  inverse  72.-1  is  only  continuous  on  S(lRr)  and  on  the  Sobolev 
spaces. 

For  polynomials  it  is  known  that 

(nC$(a.))(s,t)  = Cj!(s)cji(at),  aeSr~\  {s,t)  e Zr  (3.2) 

(see  Davison,  Griinbaum  [2])  and,  more  generally, 

(npm)(s,t)  = cji(s)pm(t),  ( s,t)eZr , (3.3) 

* 

where  the  polynomials  Pm  GlPf  (Sr~1)  are  generated  by  the  Gegenbauer  polynomials, 
i.e.  -^—[Cjifax)  = X)|m|=^  amPm(x)-  These  polynomials  Pm,  |m|  = p,  are  known  to 

constitute  a basis  for  JP^  (5r_1). 

' r 

Let  V*  :=  span{Pm  : \m\  = p).  Since  the  Gegenbauer  polynomials  CJ  can  also  be 

- * 

interpreted  as  the  reproducing  kernel  of  Mf+2(Sr+1),  the  orthogonal  projection  Fv  of 

F e C(Br)  onto  Vf  (Br)  can  be  identified  with  the  orthogonal  projection  of  F onto 
* 

Mf+2(Sr+1)  (see  Reimer  [7]  for  details).  Thus  the  theory  of  Laplace  series  can  be  used 
here  for  the  reconstruction  of  F from  its  Radon  transform. 

Let  A be  a matrix  transformation  as  introduced  in  Section  2 and  let  Fv  be  the 
orthogonal  projection  of  F onto  Vf  (Br).  Then  according  to  the  summability  theory 
of  Laplace  series  F — lim^^oo  a^F^.  Since  the  Radon  tranform  is  linear  and 
continuous  there  is  7 ZF  — lim^^oo  YX=o  anvP-Fv- 
It  can  be  shown  that  (see  Reimer  [7]) 

F„(z)  = [ (KF)  (s, t)C$  (s)d  ( tx)d(s , t),  (3.4) 

r — i Jz>- 

where 


(r  — 1)  CS  (1) 


UJr—  i * UJ-) 


2)^~ds 


2 u + r 


03, 


(3.5) 


r—  1 
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^From  this,  after  some  lengthy  calculation  using  the  adjoint  operator  of  H (which  essen- 
tially is  the  inverse  operator  of  H),  the  reconstruction  formula  follows 

F(x)  = lim  / / (7 lF)(s,t)C3  (s)Cg  (tx)dsdt.  (3.6) 

u=o  Jsr~lJ-i 

Because  of  the  identification  of  the  orthogonal  projection  of  F onto  VjJ (Br)  and  onto 
* 

J£/£+2(5r+1),  convergence  of  the  Cesaro-means  follows  for  k > r/2,  and  positivity  of 
the  operators  is  valid  for  k > r + 1.  For  the  same  reason  the  coefficients  utJU  in  the 
Newman-Shapiro  summability  method  have  to  be  calculated  for  A = f-r  1 

4 Numerical  implementation 

For  the  reconstruction  of  F formula  (3.6)  was  used.  As  soon  as  the  Radon  transform  of 
F is  known,  the  numerical  implementation  in  principle  reduces  to  a stable  evaluation 
of  the  Gegenbauer  polynomials  and  a suitable  approximation  of  the  integrals  in  (3.6). 
The  Gegenbauer  polynomials  were  evaluated  by  their  recurrence  relation  (see  Szego  [11]) 
which  is  known  to  be  numerically  very  stable.  The  coefficients  a,fU,  for  the  Cesaro-means 
and  the  Newman-Shapiro  operators  were  computed  with  the  aid  of  formula  (2.3)  and 
(2.7),  respectively.  The  factor  A„,r  was  obtained  by  (3.5).  Since  the  calculation  of  aM„  for 
the  Newman-Shapiro  operators  is  very  time  consuming  (more  than  10  hours  for  /i  > 100) 
these  coefficients  were  stored  before  the  main  computation  was  started. 

Since  the  integrand  in  (3.6)  is  a polynomial  of  degree  v + 2 with  respect  to  s (see 
(3.6)  together  with  (5.1)),  ..ds  was  approximated  by  a Gaussian-Legendre  quadrat- 
ure of  degree  p/2  + 1.  This  choice  ensures  that  for  the  evaluation  of  HF(s,t)  enough 
evaluations  with  respect  to  s are  performed  and  that  the  integral  is  evaluated  exactly 
within  numerical  precision. 

For  the  quadrature  on  Sr~ 1 first  an  interpolatory  quadrature  as  introduced  in  [6]p.l32 

was  used.  The  weights  of  such  a quadrature  formula  are  obtained  as  solutions  of  a linear 

* 

system  of  equations  GA  = e,  where  e = (1,...,1)T  e 7RV,  N = dim  7F£  (5r_1), 
A = (Ai,. . . , An)t  the  vector  of  weights  and 

G = ( XjXk ) + C^_1(xjXk))^k=1 . 

Wr-l 

The  points  were  chosen  to  be  regularly  distributed  on  latitudes  of  the  sphere. 

For  /i,  > 70  in  the  computation  of  the  weights  computational  problems  occured  be- 
cause of  a lack  of  memory.  Apart  from  this  problem,  several  weights  turned  out  to  be 
negative  which  led  to  oscillations  of  the  reconstruction.  Therefore,  this  interpolatory 
quadrature  was  substituted  by  a product-Gauss  formula  for  the  sphere  S'r"1  as  sug- 
gested by  Stroud  [10]p.  41.  The  points  and  weights  of  the  Gaussian  quadrature  were 
computed  by  the  MATLAB  program  qrule.m  which  is  available  via  internet  from  the 
Mathworks  Inc.  The  number  of  points  of  the  product  Gauss  formula  is  N — 2Mr~ 1 
where  M = p/2  + 1 is  the  number  of  points  used  in  each  direction,  i.e.  N = 2 M2  for 
r = 3. 
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All  codes  for  computation  were  written  in  MATLAB  6.  The  actual  computation 
took  place  on  a SUN  UltralO  with  256  MB  main  memory,  691  MB  virtual  memory  and 
SUN  OS  operating  system  release  5.7.  To  increase  the  computatinal  speed  all  parts  of  the 
MATLAB  code  were  written  with  as  few  for-loops  as  possible.  This  gave  an  improvement 
in  speed  of  a factor  > 500. 

5 Computational  results 

The  theoretical  results  have  been  applied  to  the  so  called  Shepp-Logan  phantom  which 
is  usually  used  as  a test  function  for  tomographic  reconstrution  algorithms.  It  is  a three 
dimensional  model  of  a human  head  consisting  of  10  ellipsoids  (see  Shepp  [9])  which  were 
shrinked  here  to  fit  into  the  unit  sphere  S2.  Figure  1 shows  a cut  at  xs  = 0.2721. 

Let  a¥\ctf\c$\  j = 1, . . . , 10,  denote  the  axes  of  the  j-th  ellipsoid,  denote  its 
density  value  and  the  diameter  of  the  ellipsoid  in  the  direction  of  t G S2 . Since 

the  Radon  transform  is  linear,  the  Radon  transform  of  the  Shepp-Logari  phantom  can 
be  calculated  to  be 

10  / (j)  ( j ) 

KF(s,t)  = ^ndb)a{3)a<23)a<3j)(s  - s[3))(s{23)  - s)  2^- 

Figure  2 shows  the  reconstruction  results  according  to  formula  (3.6)  for  Cesaro-means 
of  index  k = 4 and  for  Newman-Shapiro  operators. 

The  values  k = 1.6  and  k = 2 were  tested,  too, 
but  for  high  degrees  of  p no  convergent  beha- 
viour could  be  observed. 

For  Cesaro-means  with  k = 4 and  for  Newman- 
Shapiro  operators  Figure  2 clearly  shows  an  im- 
proving behaviour  of  the  reconstructions  for  in- 
creasing p. 

The  Newman-Shapiro  operators  show  a better 
convergence  and  for  p > 150  even  the  small 
structures  in  the  original  head  can  be  detected 
in  the  reconstruction.  It  can  be  expected  that 
for  higher  degrees  of  p this  behaviour  will  be- 
come more  evident. 


Unfortunately,  for  p > 170  the  computation  of  the  coefficients  for  the  Newman- 
Shapiro  operators  caused  some  numerical  problems  so  that  the  calculations  were  stopped 
with  p = 160.  Although  the  numerical  results  look  quite  promising,  the  drawback  in  the 
reconstruction  is  the  computational  time.  For  p = 160  the  computation  took  27.5  hours 
for  the  Radon  transform  and  31  hours  for  the  evaluation  at  the  points  x G [—  1,  l]2.  The 
evaluation  was  done  on  an  equidistant  grid  of  200  x 200  points. 


i 


0.5 


-1  -0.5  0 0.5  1 

Fig.  1.  Shepp-Logan  phantom. 
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In  principle  there  is  no  problem  to  produce  three  dimensional  reconstructions.  The  eval- 
uation points  x only  have  to  be  chosen  from  a grid  in  [— 1,1]3.  Because  of  the  time 
consuming  calculations  this  was  not  done  here,  yet. 
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